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Abstract 

The second-order factorizablc discretization of 
the compressible Euler equations developed by 
Sidilkover is extended to conservation form on gen- 
eral curvilinear body-fitted grids. The discrete 
equations are solved by symmetric collective Gauss- 
Seidel relaxation and FAS multigrid . Solutions for 
flow in a channel with Mach numbers ranging from 
0.0001 to a supercritical Mach number are shown, 
demonstrating uniform convergence rates and no 
loss of accuracy in the incompressible limit. A so- 
lution for the flow around the leading edge of a 
semi-infinite parabolic body demonstrates that the 
scheme maintains rapid convergence for a flow con- 
taining a stagnation point. 

Introduction 

Steady inviscid flow is described by the Euler equa- 
tions, which can be thought of as two subsystems. One 
subsystem corresponds to the equations governing en- 
tropy and vorticity advection. This subsystem is hyper- 
bolic in space. The other subsystem corresponds to a full 
potential operator, which is elliptic for subsonic flow and 
hyperbolic for supersonic flow. For a purely supersonic 
flow, space marching is the most efficient way of solving 
the Euler equations. For subsonic flow, the elliptic! ty of 
the full potenticil factor should be effectively handled by 
multigrid. However multigrid is not effective for advec- 
tion operators, as the coarse grid only gives part of the 
correction for certain smooth components of the error. 
Existing multigrid methods for subsonic and transonic 
flow rely on the coarse grid to smooth the entire system. 
As such, they are fundamentally limited by the ineffec- 
tiveness of the coarse grid in correcting the part of the 
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error corresponding to advection factor. This same diffi- 
culty is true for high Reynolds-number viscous flows as 
well. 

As has been pointed out by Brandt^, to obtain ideal 
multigrid convergence rates for subsonic, inviscid flows, 
the discretization must distinguish those parts of the dif- 
ferential operator which correspond to advection, and 
those which correspond to elliptic behavior. The advec- 
tion terms are treated efficiently by marching and the el- 
liptic terms are rapidly solved with a multigrid iteration. 
Brandt * argues that by splitting the system into its ellip- 
tic and advection parts, the convergence rate of the full 
system ought to be equal to the slowest of the two sub- 
systems. Using this approach, Brandt and Yavneh have 
demonstrated textbook multigrid for the incompressible 
Navier-Stokes equations^ . Their results are for a simple 
geometry and a Cartesian grid, using a staggered-grid 
discretization of the equations. In a closely related ap- 
proacli, Ta’asan presented a fast multigrid solver for the 
compressible Euler equations. This method is based on 
a set of “canonical variables” which express the steady 
Euler equations in terms of an elliptic and a hyperbolic 
partition^. In Reference 3 it is shown that ideal multi- 
grid efficiency can be achieved for the compressible Euler 
equations for two-dimensional subsonic flow using body- 
fitted grids. 

In an earlier work^ the authors presented a multi- 
grid scheme for the steady, incompressible Euler equa- 
tions based on a pressure Poisson discretization which 
distinguishes the advection operator from the elliptic 
part of the system of equations. Both structured grid 
and unstructured grid flow solvers were written using 
the discretization. The results in Ref. 5 demonstrated 
that the scheme can achieve ideal multigrid convergence 
rates for internal flows. The scheme has been extended 
to three dimensions by Sanchez , who has also demon- 
strated ideal multi grid convergence rates for internal 
flows. This approach has the advantage of non-staggered 
grids a collocated, vertex-based discretization of the 
equations is used, simplifying the restriction and prolon- 
gation operations and allowing the use of simple point 
collective Gauss- Seidel relaxation. Although the pres- 
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sure Poisson approach may be extended to compressible 
flows, it is not conservative and is unsuited for super- 
critical flows with shocks. Furthermore, the extension 
to viscous flows is limited to the incompressible Navier- 
Stokes equations. 

Recently, Sidilkover 0 has devised a discretization of 
the compressible flow equations that overcomes these 
limitations. This discretization may be applied to the 
Euler equations in conservative form, using the multidi- 
mensional upwind scheme of Sidilkover^’ In Ref. 10 
it is shown that this approach should lead to a scheme 
that is factorizable , i.e., the scheme distinguishes be- 
tween those parts of the operator that represent advec- 
tion, and that part of the operator that represents po- 
tential flow. In Ref. 8, such a factorizable scheme is 
constructed for Cartesian grids. Because the discretiza- 
tion is stable for Gauss-Seidel relaxation, it the conver- 
gence rate does not depend upon a Courant-Friedrich- 
Lewy number restriction, unlike standard discretizations 
of the Euler equations which must use time marching. 
For this reason, the same convergence rate is obtained 
for subsonic Mach numbers all the way to the incom- 
pressible limit. Tlie scheme may be written in the form 
of a central- difference part plus an artificial viscosity. 
As such, it is very similar in formulation to standard 
upwind, finite- volume discretizations for the Euler and 
Reynolds- Averaged Navi er- Stokes equations, and can be 

written as a conventional upwind scheme with a rnodi- 
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fled artificial dissipation. Sidilkover shows that for the 
discrete scheme to preserve factorizability the dissipation 
terms must be discretized in a specific way. In addition, 
he shows that the artificial viscosity may be rescaled with 
the Mach number such that the factorizability, and thus 
the accuracy, as well as the h-ellipticity of the operator 
is preserved in the incompressible limit. 

In the present work, the generalization of Sidilkover’s 
factorizable scheme to curvilinear, body-fitted coordi- 
nates is presented. First, the governing equations are 
presented, including the form of the multidimensional 
upwind artificial dissipation. It is shown how the dissi- 
pation terms must be discretized to maintain factoriz- 
ability. A discussion of the point collective Gauss-Seidel 
relaxation is presented next . Solutions for flow in a chan- 
nel, with inlet Mach numbers ranging from 0.0001 to 
a supercritical Mach number are shown, demonstrating 
the accuracy and convergence rates of the scheme. An 
additional computation for the flow around the leading 
edge of a semi-infinite parabola demonstrates that, the 
current scheme does not suffer from the convergence dif- 
ficulties near a stagnation point that were observed for 
the authors’ previous scheme 

Mathematical Formulation 

The artificial dissipation of the factorizable scheme 
can be described by first presenting the modified equa- 
tion, or first differential approximation (FDA), of the 
discrete scheme. This is the differential equation which 
is found by expanding the difference equat ion in a Taylor 


series about each grid vertex and considering the lead- 
ing terms of the truncation error. These terms are the 
artificial dissipation of the scheme. 

The starting point for the scheme is the two- 
dimensional Euler equations in non-conservation form. 
Let p be the density, u = tu + jv be the velocity, and p 
be the pressure. The entropy s is defined as 


•■(£)(£)" 

(i) 

where po and po are a reference pressure and density, 
respectively, and 7 is the ratio of the specific heats. 
Then the Euler equations may be written in the vari- 
ables (s, w,r,p): 

?1Vs = 0, 

(2a) 

u • v a + i vp = o, 

p 

(2b) 

pc 2 V*u -1- u-Vp = 0. 

(2c) 


The factorizability of the scheme depends on the form 
of the artificial dissipation added to this system of equa- 
tions. 

The entropy is only weakly coupled to the momen- 
tum and the pressure equations through the equation 
of state (1). In fact, the entropy equation corresponds 
to one of the advection factors of the Euler equations. 
Therefore, it may be discretized independently of the 
momentum and pressure equations in any appropriate 
way without affecting the factorizability of the scheme. 
The advection operator u V uses simple upwind differ- 
encing in Eq. (2a). Let (£,17) be a general curvilinear 
coordinate system and define the contravariant compo- 
nents of the velocity, ( U , V) by the transformation 

G: :;)(?)=(“)- 131 

In this coordinate system tTVff = Ud$ + Vd v * The equa- 
tions are discretized on a uniform grid in (£ T f?) space, 
with a grid spacing A£ = Ap = 1. The FDA of the 
first-order upwind difference approximation to uV is 

q = «V-i|£7|^-i|V-|c^ = 0 , (4) 

and the entropy equation is discretized as 

qs = 0. (5) 

A second-order upwind discretization of the advection 
operator has also been used in Eq. ( 5 ). However, for 
the particular cases shown below the use of a second- 
order advection operator has an insignificant affect on 
the computed results. This point is discussed in the Re- 
sults section. 

The dissipation for the momentum and pressure equa- 
tions. Eqs. (2b) and ( 2 c), is the multidimensional upwind 
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dissipation of Sidilkover^. In vector notation, this dissi- 
pation is written as 

u-V«7+ -Vp- —V (pr 2 V u + u Vp) = 0, (6) 

p 2pc ' 

pc 2 Vu + t7Vp — pc^-V- (u-Vu + =0, (7) 

where c is the speed of sound, <r m and cr p are scaling 
coefficients, and C is a length scale proportional to the 
grid spacing. Note that the dissipation of the momen- 
tum equation is the gradient of the pressure equation 
residual, and the dissipation of the pressure equation is 
the divergence of the momentum equation residual. It. 
is this property of the artificial dissipation that leads to 
a factorizahle scheme. Note that the curl of Eq. (6) is 
identical to the curl of Eq. (2b), i.e., vorticity equation of 
the governing Euler equations is unaffected by the arti- 
ficial dissipation. The vorticity equation corresponds to 
the second advection factor of the Euler equations. Sim- 
ilarly, because the dissipation of Eq. ( 7 ) is the divergence 
of the momentum equation, it is affected only by the ir- 
rotational part of the velocity field but not the solenoidal 
part. The irrotational part may be written as the gradi- 
ent of a potential. As the pressure equation Eq. ( 2 c) is 
a form of the continuity equation, it corresponds to the 
full potential factor, and the FDA described by Eq. ( 7 ) 
preserves this property. 

If the advection operator, pressure gradient, and di- 
vergence terms in Eqs. (6) and ( 7 ) are discretized us- 
ing central differences, the scheme is second-order accu- 
rate and factorizable. However, such a scheme is not. 
h-elliptic. This lack of h-ellipticity is a result of the 
central difference approximation to the advection term 
in the momentum equation, Eq. (6). Sidilkover shows 
that this advection operator corresponds to vorticity 

o 

advection 0 . Replacing this operator by the first-order 
upwind approximation in Eq. ( 4 ) restores h-ellipticity 
and the scheme remains factorizable, but it now becomes 
only first-order accurate. Note that the pressure advec- 
tion term in Eq. ( 7 ) continues to be approximated by 
second-order central differences. 

To obtain second-order accuracy while maintaining 
h-ellipticity, appropriate antidissipative term? must be 
added to the momentum equation in such a w r ay that fac- 
torizability is preserved. The form of these terms is de- 
pendent upon the computational coordinates (£,rj), and 
they must be written in terms of the contravariant and 
covariant components of the velocity vector. The covari- 
ant components are related to the physical components 
by the transformation 



Writing the advection operator in Eq. (6) in terms of 
the covariant velocity components and ignoring terms 
containing the higher-order geometric derivatives, the 
scheme may be upgraded to nearly second-order by 


adding the following corrections: 

u VU = qU + |f'| d\u + i \V\ d ( c)„v, ( 9 a) 

u-VV = qV + i \U\ dtd„U + i |V| t) 2 V. ( 9 b) 

The cross-derivative terms in Eq. ( 9 ) are necessary for 
second-order accuracy and to retain factorizability at the 
same time. Using Eq. (8) to find U and V, the antidis- 
sipative terms in Eq. ( 9 a) and ( 9 b) are evaluated, and 
then rotated back into the physical components. This 
gives a discretization of the form 

tFV« = q« + VD (10) 

w T here 

D= l -(\ll\d ( U + \V\d v V). ( 11 ) 

The above expression is the generalization to curvilinear 
coordinates of the one given in Ref. 8. 

To gain more insight about the nature of the correc- 
tion terms, the FDA of Eq. (10) may be rewritten as 

qt 7 +VD = u-Vu + J - (r ( \V \ rV - e" |P| d<«) 

where = d x v — d y u is the vorticity, J = 
is the Jacobian of the coordinate transformation, and 
(e^c 75 ) are the contravariant basis vectors in the (£,77) 
coordinate directions. Note that, if the flow is irrota- 
tional, then the first-order truncation error terms vanish 
identically and the approximation becomes second-order 
accurate. 

It was pointed out above that factorizability depends 
upon the dissipation of the momentum equation being 
written as the gradient of the pressure equation, and the 
dissipation of the pressure equation being written as the 
divergence of the momentum equation. Likewise, the 
second- order correction to the advection terms must be 
in the form of a gradient so that the vort icity equation 
is is unaffected. With the definition in Eq. (10), taking 
the curl of Eq. (6) yields qu) = 0. Thus the advection 
operator acting on the vorticity is unchanged from the 
first-order scheme. 

To get full second-order accuracy it is necessary to 
use second-order accurate discretization of the advection 
operator in the momentum equation. The FDA of an h- 
elliptic, fully second-order upwind advection operator 

qHo = Ud$ + Vd r] 

- i (|C| d\ + 2 sgnP Vd ( d v + ( 12 a) 

when |t/| > |U| and 

quo — Udz + Vdr, 

- \ ' 9? + 2 S6 nr Ud(d n + \v\ dl , j (12b) 

* = «y» 7 - jx,„ J ?’ 1 = -iy^ + jx 4 
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when |C?| < | V | . With this advection operator in the 
momentum equation, the quantity D in Eq. (11) must 
be modified in order to preserve the factorizability of the 
scheme. This modified D is given by 


Duo = D + -sgnZJ V {drfl + c^V) 




d„V ( 13 a) 


when | U | > |l' r | and 

Dho = D + isgnF U (d v U + 3 { V) 


\U ( 13 b) 


when | f 7 1 < 1 1 7 1 . Taking the curl of the momentum 
equation now gives qHO w = 0, i. e., the vorticity equa- 
tion is now second-order accurate. 

To see how the multidimensional upwind dissipation of 
Eqs. (6), ( 7 ) and ( 9 ) is related to the standard first-order 
upwind difference scheme, consider a uniform Cartesian 
grid with equal spacing in the x and y directions. If the 
scaling coefficients <r m and <J p are taken to be one, i is 
the grid spacing, and the cross- derivative terms in the 
dissipation terms of Eqs. (6) and ( 7 ) are ignored, the 
standard first-order upwind scheme is obtained. The 
first-order upwind scheme on a nonuniform grid further 
replaces fV by i x dx + where and ( y are the grid 

spacings in the two coordinate directions. The multidi- 
mensional upwind scheme uses a single length scale, and 
retains the cross-derivative terms. Consider the advec- 
tion operator of Eq. ( 9 ) on a Cartesian gnd. Dropping 
the cross-derivative terms yields the standard first-order 
upwind discretization. 

Discretization 

The factorizability of the FDA is a necessary condi- 
tion for the factorizability of the difference scheme, but 
it is not sufficient. For the difference scheme to be fac- 
torizable, the difference operators must commute in the 

o 

same way as the differential operat ors . Introducing the 
difference approximations to the partial derivative oper- 
ators, 


d £ — $4 -+■ * * • , df} — -f- ■ * * j 

= c?| + * ■ • , ^ 4 — » 

= »! + ••• i 

the following conditions must hold 

44 , = 

d^di = 


The following difference operators satisfy this condition: 
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To write the complete discrete scheme, the subscript h 
is used to denote a standard difference to the corre- 
sponding operator, and the addition of the overbar ( ) 
denotes the “wide” differences given above. The second- 
difference expressions may be expressed in flux form by 
taking a six-point difference centered on an edge between 
two vertices, and then taking a two-point difference of 
those expressions centered on a vertex. The subscripts e 
and v are used to denote difference operators centered 
on an edge or a vertex, respectively. The fully discrete 
scheme is then written as 


= 0, 


( 14 a) 


q £o3 + V„D£ 0 + iv h p 

P 

— fpc 2 V£-u + u-V^p] = 0, (14b) 

2 pc V / 


pc 2 V h -u + nV h p 

- pc^-Vl- (t7V*« + ^ (v* + V*) p) = 0. (14c) 


The derivatives in that part of Dho given in Eq. ( 11 ) 
are discretized using the wide, six-point difference stencil 
for d^U and d*V. The additional terms of Dho and 
the corresponding dissipation terms for the second-order 
advection operator qHO must be discretized in a very 
precise way, and depend upon the flow direction and 
the relative magnitudes of U and V. The details are 
described in Ref. 12 for uniform Cartesian grids. The 
stencils for the general curvilinear grids are presented in 
the appendix to this paper. 

The scaling coefficients are 


a m — max (M, Af c ), 


1 

max (A/, M c ) 


( 15 ) 


where M is the local Mach number, and A/ c is a cutoff 
Mach number to prevent division by zero. The cutofT is 
chosen to be 0 (h ), and essentially becomes active near 
stagnation points. The purpose of the rescaling of the 
pressure equation dissipation, is to prevent the ellip- 
tic factor in the discrete equations from becoming the 
skewed Laplacian operator in the incompressible limit. 
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Currently, we take M c = i' | A A/ j , where z> = 5 
and A A/ is the two-point, difference in Mach number 
on an edge. For the channel flow cases shown below, A/ c 
never becomes larger than M. For the leading-edge flow, 
the cutoff does become active at the stagnation point. 
When M c > M, it is necessary to add additional dissi- 
pation to the advection operator qjj 0 . The form of this 
dissipation is a five-point, pseudo- Laplacian, 

d . P = \ max (0, M c - M) j (d\ + dl) , (10) 

where J is the Jacobian of the coordinate transforma- 
tion. This operator is added to both the entropy and the 
momentum equation, and is cast in flux form in order to 
maintain conservation. No attempt has been made to 
optimize either the coefficient v or the form of the oper- 
ator d sp . 

As written, the Eq. (14) is valid only for subsonic 
flows. This is because the pressure differences in the 
artificial dissipation terms are not fully upw'inded in a 
supersonic zone. A simple modification to Eq. (14) can 
be made by rescaling the gradients of the pressure when 
the flow becomes sonic. Introducing the parameter k, 
defined as 


k = max(l, A/ 2 ), (17) 

the final form of the scheme is 

q h 5 - d sp s = 0, (18a) 

qfiow - d 3 pu + V„Dho + ~^ h P 

- ^ (pc 2 ?*i? + ifl-Vjp) = 0, (18b) 


operators in Eq. (14). The artificial dissipation in 
Eqs. (14b) and (14c) can be rewritten as 


Su = — Vt-Dno + d S p(i 


Sp = 




+ vj 



[¥l 

cVg-u-h — tT^p^l 

. (20) 

2 ' 

V, PCK J J 

-» r-? h 

peu ■ V 

S+f (vf + «)»)] 

• (21) 


This is a conservative form of the dissipation cast in 
terms of the primitive variables. The terms inside the 
square bracket in Eqs. (20) and (21) are now interpreted 
as dissipative fluxes. These terms are discretized on cell 
faces according to the appropriate wide or narrow differ- 
ences. The first-order upwind advection operator q ; ' and 
the antidissipative terms D h (U,V) in Eq. (20) may be 
recast in conservation form and evaluated on cell faces. 
The length scale ^ is evaluated on the cell face by tak- 
ing min(^, 4?), the shortest length in each grid direction. 
The scaling coefficients <r m and <r p are evaluated using 
the Mach number on the face. 

In addition to the dissipation terms, examination of 
Eqs. (14b) and (14c) sho^vs that the pressure gradi- 
ent terms are discretized using the wide stencils. The 
central-difference part of the conservation equations (19) 
must be corrected to account for these wide differences. 
This is done by adding a term to the dissipative fluxes 
for the momentum and pressure equations as follows: 


6u<- 53- - (V - V h ) p, 

(22) 

p\ / 

6p Sp - u- (v h - p. 

(23) 


Once the dissipative fluxes (As, Att, Ac, Ap) have been eval- 
uated, they are converted to the conservation variables 
by the transformation 


pc 2 V h u -j- u-V h p 

- (tT-VeU + ^ (Ve + V?) p) = 0. (18c) 


( &P \ 


-P/s 

0 

0 

!/c 2 \ , 

( Ss\ 

Spu 


-pujs 

p 

0 

u/c 


Su 

Apt? 


-pv/s 

0 

P 

t'/c 2 


Sv 



\-p(u 2 + v 2 ) /2s 

pu 

pv 

h/c 2 ) 


VpJ 


(24) 


Because the difference equations (14a), (14b) 

and (14c), can be written as a central difference part 
plus a dissipation, it is straightforward to obtain a con- 
servative discretization. The conservation form of the 
Euler equations are discretized using a central- difference 
finite- volume approximation, 


v h (p »0 = o, 

(19a) 

V h (puit) 4- V h p = 0, 

(19b) 

V h [(p f +p)t7] = 0 

(19c) 


where e = u*u/2 + p/(p(l — 1)) is the total energy. 
The dissipative fluxes on each cell face are computed 
in terms of (s,u,u,p) using the appropriate difference 


w'here h = e + p/p is the total enthalpy. 

Solution Procedure 

The solution of the discrete equations is performed us- 
ing a symmetric point collective Gauss-Seidel iteration, 
which has been found to be a very effective smoother. 
The grid points are ordered such that, the forward Gauss- 
Seidel sweep is in the streamwdse direction. This means 
that the entropy equation is marched in space during the 
forward sweep. The reverse Gauss-Seidel sw r eep is under- 
relaxed with a factor of 0.5 for stability. The residuals 
of the conservation equations (r p , r pu , r pV} r pe ) are com- 
puted at each vertex using the discretization presented 
in the previous section. These residuals are then trans- 
formed to the residuals of the primitive variables by the 
inverse of the transformation in Eq. (24). 
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Figure 1. Channel geometry. 


At each point, the update to the solution is given by 


M 


( a-A 

Am 

A,- 

Vp)„ 



(25) 


where M is the matrix of coefficients of the primitive 
variables at vertex The coefficients are found by 

collecting the contributions to vertex (i, j) from the dissi- 
pation terms on the four surrounding faces. Because the 
entropy equation decouples from the rest of the system, 
this is a block diagonal matrix where the upper block is 
the upper left-hand entry, and the lower block is a 3 X 3 
matrix of coefficients multiplying v tJ . and p l;J . This 
matrix is easily inverted. 

The relaxation is accelerated using a standard Full- 
Approximation Scheme (FAS) multigrid. A sequence of 
grids (7 a' j G/v _ i , . . . ,Cj q is used, where (7 a is the finest 
and Go the coarsest. Let Lfc-i be the coarse grid oper- 
ator, u be the vector of the conservation variables, 
be the fine- to- coarse grid restriction operator, and 7* _1 
be the coarse- to- fine grid prolongat ion operator. If Uk is 
the current solution on grid k> the residual on this gnd 
is r* = ffc — LfrU*. This is the residual of the conserva- 
tion equations, not the primitive equations. This leads 
to the coarse-grid equation 

It-iu*-, = r*_i = /£_,r* + Lfc-i (/*_,U*) . ( 26 ) 


After solving the coarse-grid equation for Uk-u the fine- 
grid solution is corrected by 

ur" «- u k + /£- 1 (u*_, - /£_ ,fu) . (27) 


Equation (26) is solved by applying the same relaxation 
procedure that is used to solve the fine-grid equation. 
Multigrid is applied recursively to the coarse-grid equa- 
tion. On the coarsest grid, many relaxation sweeps are 
performed to insure that the equation is solved com- 
pletely. A conventional V cycle or IF -cycle is used. 


Results 

Solutions are shown for compressible flow in a two- 
dimensional channel, the geometry of which is shown in 
Fig. 1 . The shape of the low^er wall between 0 < x < 1 
is t/(:r) = rsin 2 7nr. For the computations shown here, 
the thickness ratio r is 0.05. The grid spacing is uniform 
in the ar-direction, and the coordinates in the y - direction 
are found using a simple shearing transformation. 



Figure 2. Mach number contours for a M = 0.5 
inlet., contour increment AM = 0.01, for a 385 x 129 
grid. 


At the inflow* boundary the entropy, total enthalpy 
and flow* angle are specified, and the pressure is extrapo- 
lated from the interior. The outflow* boundary condition 
is a specified pressure and s, u and v are extrapolated 
from the interior. At the upper and lower w r alls of the 
channel, the flow* tangency condition u n — 0 is enforced 
by setting the residual of the momentum equation nor- 
mal to the w*all to zero. Because of the wide stencils the 
dissipation terms on the wall are evaluated using a row* 
of ghost vertices. The pressure is extrapolated to those 
vertices using the normal momentum equation. The en- 
tropy and total enthalpy at the ghost vertices are set to 
the values on the w'all, and the normal component of the 
velocity is reflected from the interior. 

Solutions arc obtained using a FMG cycle to initial- 
ize the solution. A solution is computed on the coarsest 
grid and prolonged to the next grid to obtain a starting 
solution for that grid. This procedure is continued re- 
cursively until the finest grid is reached. On each grid, 
a Vfl, 1) multigrid cycle is used to solve the system. 
Five multigrid cycles are run on each of the coarse grids, 
and fifteen cycles are run on the finest grid. After each 
symmetric Gauss-Seidel relaxation sw r eep, an additional 
three streamwise sweeps are done on the w*all and its 
neighboring rowr of vertices. 

Mach contours are seen in Fig. 2 on a 385 x 129 fine 
grid for an inlet Mach number of 0.5. The solution 
is symmetric except for a glitch at the outflow bound- 
ary, which is caused by specifying uniform pressure at 
the outflow’ boundary. Convergence rates are shown in 
Fig. 3. A total of 7 grids is used, the coarsest being 7x3. 
The residual levels are renormalized after each prolonga- 
tion to the next finest grid, and on the coarsest gnd 
several relaxation sw 7 eeps are performed, giving essen- 
tially a direct solve on that grid. The convergence rate 
on the finest grid is initially 0.25 per V'(l, 1) cycle, and 
the asymptotic rate is 0.38 per cycle. The residual re- 
duction over 5 cycles is seen to be the same on each grid, 
showing that the convergence rate is O(n). The drag on 
the lower w*all and the rms error in the entropy over the 
computational domain are showm in Fig. 4. On the finest 
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Figure 3. Convergence rates for a M = 0.5 inlet flow 
using a FMG cycle with a 385 x 129 fine grid. 


grids both are converged after essentially two multigrid 
cycles. 

The drag and entropy appear to be exhibiting bet- 
ter than first-order accuracy but not quite second-order 
accuracy. On the finer grids, the error decreases by ap- 
proximately one-third with the halving of the grid spac- 
ing. The entropy error decreases by about 0.28 on the 
finest grids, while the drag is reduced by about 0.35. 
These values are somewhat sensitive to the choice of the 
wall boundary condition. Curiously, when simple reflec- 
tion boundary conditions are used, the drag converges 
better than second order, dropping by a factor of six 
to eight with each grid refinement. Neither the drag 
nor the entropy is very sensitive to whether a first-order 
or a second-order advection operator is used in either 
Eq. (18a) or Eq. (18b). This is because the flow is a 
potential flow, and the factorizability of the scheme de- 
couples the entropy and vorticit.y advection from the po- 
tential factor of the operator. Given that s and jj are 
identically zero analytically, first-order advection is per- 
fectly adequate. 

To demonstrate that the factorizable scheme requires 
no special preconditioning to handle very low Mach num- 
bers, a solution for an inlet Mach number of 10“ 4 is 
shown in Figs. 5, 6 and 7. Because the ratio of the 
dynamic pressure to the static pressure is 0(A/ 2 ), the 
relative roundoff error is much larger for this case than 
for the higher Mach number cases, so the residuals bot- 
tom out at a higher level. The initial convergence rate 
is 0.24 per cycle, very close to that for the Mach num- 
ber 0.5 case. As with the higher Mach number case, 
the convergence rate is G(n ). The solution is seen to 
be very symmetrical at the bump, with some boundary 
layer behavior at the inflow and the same behavior at 
the outflow as for the the higher Mach number. The 



Figure 4. Convergence of the Ln entropy error and 
drag on the lower wall for a M = 0.5 inlet flow using 
a FMG cycle with a 385 x 129 fine grid. 

drag and entropy errors converge within two cycles on 
the finest grids. 

Mach contours, residual convergence and the drag and 
entropy are shown in Figs. 8, 9 and 10 for an inlet Mach 
number of 0.73, which corresponds to a supercritical inlet 
Mach number. The peak Mach number before the shock 
is approximately 1.2. The rate of residual reduction is 
initially about 0.53 per cycle on the 385 x 129 fine grid, 
and the asymptotic rate is 0.68. The drag and entropy 
take longer to converge than for the subcritical cases, 
and in fact they have not converged oil any of the coarse 
grids before the solution is interpolated to the next finer 
grid. On the finest grid, the drag and entropy converge 
in about eight cycles. 

All the solutions shown so far do not have a stagna- 
tion point. In earlier work, the pressure Poisson scheme 
of Roberts, Sidilkover and Swanson 0, exhibited a de- 
terioration in the convergence rate for flows containing 
a stagnation point ^ ^ . The current scheme does not suf- 
fer from this difficultly. A solution was obtained for the 
flow about the leading edge of a semi-infinite parabolic 
body, shown in Fig. 11. For this case, Dirichlet boundary 
conditions were used at the far-field, where the velocity 
field was taken from the incompressible solution given by 
a conformal mapping around the body. The freestream 
entropy was specified, and a uniform total enthalpy was 
computed by taking the freestream Mach number equal 
to 0.1. A fine grid of 129 x 129 grid points was computed 
from the conformal mapping, and 6 grid levels were used. 
In the symmetric Gauss-Seidel relaxation, it was neces- 
sary to underrelax the stream wise sweep with a factor 
of 0.9 for stability. The wall boundary conditions were 
simple reflection of the velocity, entropy and pressure. 

Pressure coefficient contours are shown Fig. 12, and 
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Figure 5. Mach number contours for a M = 10~ 4 
inlet, contour increment AM = 2x 10~ 6 , for a 385 x 
129 grid. 



Figure 6. Convergence rates for a M =10 4 inlet 
flow using a FMG cycle with a 385 x 129 fine grid. 


are seen to be symmetric about the axis of the parabola. 
Convergence rates are shown in Fig. 13. Unlike the chan- 
nel flow, the convergence rates are not quite O(n), be- 
coming slightly slower on finer grids. This is because 
of the underrelaxation of the streamwise Gauss-Seidel 
sweep, which prevents the advection factor from being 
solved in that sweep. The asymptotic convergence rate 
on the finest grid is approximately 0.50 per cycle. The 
entropy error, Fig. 14, is seen to converge as well as for 
the channel flow, with a reduction of aronnd one-third 
with each grid doubling. 

Conclusions 

A factorizable discretization of the compressible Eu- 
ler equations on general curvilinear body-fitted grids has 
been presented. The discretization is based on the mul- 
tidimensional upwind scheme of Sidilkover, and can be 
written in a form that is closely related to conventional 
upwind-differenced finite- volume schemes. Unlike con- 
ventional schemes, this discretization lends itself to an 
extremely efficient multigrid solution procedure. The 
factorizability of the scheme has the property that the 
correct incompressible limit of the compressible equa- 
tions is obtained without the need for any special pre- 
conditioning. Because the scheme is stable for Gauss- 
Seidel relaxation, fast convergence may be obtained for 
essentially incompressible flow as well as high subsonic 
Mach numbers. This discretization unifies compressible 
and incompressible flow algorithms. 
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Appendix: Discretization of D ho 

The discretization of the derivatives in the expression 
for Dho, given in Eqs. (13) must be done in a very precise 
w 7 as so as to preserve the factorizability of the discrete 
scheme. In Ref. 12, Sidilkover presents these expres- 
sions for a uniform Cartesian grid. The particular case 
of |t7| > [ V | is presented here. In this case, Dho is given 
by Eq. (13a), repeated here for convenience. 


Third, take U <0, V > 0. 

C>*V = \ (Vi+2,2 + V.+2,i-l - V.+ i,, - V.+ .,,-. ) . 
^V=i(V 1+2 , J+l -V.+2,,-,). 

Fourth, take U < 0, V < 0. 

— “ (^i+lj+2 — St + ij) , 

= - (Vt + 2 J +1 + V\-42 fJ - V.+ lJ-fl - Pi+l,j) , 

\ (Vi+2,j+l - Vi+2,2-0 . 

Now 7 consider the rj-cel] edge, (i,j A 1/2). First w 7 e 
take U > 0, V > 0, in which case the following difference 
formulas are used. 

d*U = ~ (ft,,; - ft,,,., + ft,_,,, - ft, , 

= \ (V .,7 - V’,-2,2) , 

d!;v = i (Vi-1,2+1 - Vi — 1 ,2 + Vi— 2 , 2 + 1 - Vi— 2 , 2 ) . 
Second, take U > 0, V < 0. 


Duo = D+ isgnft V' (c?„ft + d ( V) 

)ft>V 

The term D has already been discussed in the Discretiza- 
tion section, and what remains is d v U, d(V and d^V 
above. Each of these terms is discretized differently on 
a £-ce!l edge, (i ± 1/2, j), and an rj-cell edge, (i, j± 1/2), 
w r here i and j are the indices of the vertices in the £ and r/ 
directions, respectively. 

Consider the £-cell edge, (i+1/2, jf). First w r e take U > 
0, V > 0, in w r hich case the following difference formulas 
are used. 



dqU = — (UiJ+2 +^j-1 3 j + 2 —Ut- I,j+l) , 

9 { V = ^ (Vi, J +1 - V._ 2 , 2 +i) , 

SjV = \ (V,— ,.2+, - Vi-, ,2 + Vi— 2,2 + 1 - Vi— 2,2) . 
Third, take U < 0, V > 0. 

a*ft = \ (ft.+ 1,2 - ft.+i,j-i +ft.,2 2-1) , 

d?V = \ (Vi+2,2 - Vi ,2 ) , 

d*V = \ (Vi+ 2 , 2+1 - Vi+ 2,2 + V, + 2 , 2 + 1 - V+,, 2 ) . 

Fourth, take V < 0, V < 0. 


d*U = i (ft, ,2 - ft, ,2-2) , 

3? V = i ( V , ,2 + V, , 2—1 - V, — , , 2 ' - Vi— ,,j_ 1 ) , 
<V=t(V,_,,2+l -V,_,,2-l). 


4V=t(V, + 2,2+l -V..2+,), 

SjV = i (V. + 2.2+. - V,+2,2 + Vi+2,2+1 - Vi+ 1 , 2 ) - 


Second, take V > 0, U < 0. 

d*U = I (ft, , 2+2 - ft,, 2) , 

9 ? V = i (V,,2+. + Vi, 2 - V, — , ,2 + 1 - V, — 1,2) . 

K V = - (Vi — , , 2 + ! - V,: — ] ,2 — 1 ) ■ 


The above differences are also the ones used in the 
definition of the second-order advection operator qHO- 
Analogous expressions may be developed for Eq. (13b), 
corresponding to [ V | > \0\. 
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Figure 7. Convergence of the 7 2 entropy error and 
drag on the lower wall for a M = 10 -4 inlet flow 
using a FMG cycle with a 385 x 129 fine grid. 



Figure 9. Convergence rates for a M = 0.73 inlet 
flow using a FMG cycle with a 385 x 129 fine grid. 



Figure 8. Mach number contours for a M = 0.73 
inlet, contour increment AM = 0.025, for a 385 x 129 
grid. 



Figure 10. Convergence of the L<t entropy error and 
drag on the lower wall for a M — 0.73 inlet flow 
using a FMG cycle with a 385 x 129 fine grid. 
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V(1,1) cycles 

Figure 13. Convergence rates for a M = 0.1 leading- 
edge using a FMG cycle with a 129 x 129 fine grid. 

Figure 11. Semi-infinite parabola, 33 x 33 grid. 
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Figure 14. Convergence of the Lo entropy error for a 
Figure 12. Pressure coefficient contours around a M = 0.1 leading-edge flow using a FMG cycle with 

M = 0.1 leading-edge, increment AC P = 0.05, for a 129 x 129 fine grid, 

a 129 x 129 grid. 
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